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A planar crack generically segments into an array of “daughter cracks” shaped as tilted facets when 
loaded with both a tensile stress normal to the crack plane (mode I) and a shear stress parallel to 
the crack front (mode III). We investigate facet propagation and coarsening using in-situ microscopy 
observations of fracture surfaces at different stages of quasi-static mixed-mode crack propagation 
and phase-field simulations. The results demonstrate that the bifurcation from propagating planar 
to segmented crack front is strongly subcritical, reconciling previous theoretical predictions of linear 
stability analysis with experimental observations. They further show that facet coarsening is a self¬ 
similar process driven by a spatial period-doubling instability of facet arrays with a growth rate 
dependent on mode mixity. Those results have important implications for understanding the failure 
of a wide range of materials. 

PACS numbers: 62.20.Mk, 46.50.+a, 46.15.x 


Crack propagation is a main mode of materials failure. 
Understanding and controlling this complex phenomenon 
continues to pose both fundamental and practical chal¬ 
lenges. While quasi-static planar crack growth with a 
tensile stress normal to the fracture plane (mode I) is 
well-understood, geometrically much more intricate crack 
patterns can form in varied conditions [1]. A few exam¬ 
ples include thermal or drying stresses that can cause 
cracks to oscillate and branch [2, 3], or re-organize into 
complex three-dimensional patterns [4-6], nonlinear elas¬ 
tic effects that can induce crack front instabilities even 
in mode I [7], or the superposition of mode I and a shear 
stress parallel to the crack front (mode III). This mixed¬ 
mode I-j-III fracture is observed in a wide range of en¬ 
gineering and geological materials to produce arrays of 
daughter cracks, which are shaped as tilted facets and 
form by a geometrically complex crack front segmenta¬ 
tion process [8-23]. 

Recent theoretical progress has been made to charac¬ 
terize the crack-front instability leading to segmentation 
[24, 25] and to describe the propagation of daughter-crack 
arrays [26]. However, theory and experiments have not 
produced a consistent picture. Griffith’s energetic cri¬ 
terion [27] predicts that planar crack growth is possible 
when the elastic energy release rate 

^ ^ ((1 “ + ^ni) 5 ( 1 ) 

exceeds a critical material-dependent threshold Gc, 
where Kj and Km are the mode I and mode III stress 
intensity factors (SIF), respectively, which characterize 
stress divergences near the crack front, // is the shear 
modulus and u is Poisson’s ratio. Phase-field simulations 
of brittle mixed-mode I-flll fracture have revealed that 
planar growth is linearly unstable against helical defor¬ 


mations of the crack front, which couple in-plane and 
out-of-plane perturbations and develop nonlinearly into 
facets [24]. A subsequent linear stability analysis in the 
framework of linear elastic fracture mechanics (LEFM) 
[25] has predicted that this helical instability should oc¬ 
cur when Km /Kj exceeds a threshold 

Kiii \ ^ / (l-i^)(2-30 
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which only depends on Poisson’s ratio. However, 
crack front segmentation is experimentally observed for 
KmjKi values much smaller than this threshold [8, 23], 
or even vanishingly small [22]. This apparent disagree¬ 
ment between linear stability analysis and experiment 
raises the question of whether LEFM and phase-field 
modeling are adequate theories to describe crack propa¬ 
gation in mixed-mode I+HI brittle fracture. Also poorly 
understood is “facet coarsening”, the progressive increase 
of facet width and spacing with propagation length from 
the parent crack. Phase-field modeling [24] and experi¬ 
ments [28] suggest that coarsening may be a self-similar 
process, but its precise mechanism and dependence on 
mode mixity are not well understood. 

In this letter, we investigate both facet propagation 
and coarsening by mixed-mode I+HI fracture experi¬ 
ments that allow us to visualize in-situ complex crack 
morphologies during quasi-static propagation, thereby 
providing much more detailed geometrical information 
on crack front evolution than conventional post-mortem 
fractography. Moreover, we carry out phase-field simula¬ 
tions of those experiments that allow us to relate exper¬ 
imental observations to LEFM theory. The results help 
resolve the puzzling discrepancy between linear stability 
analysis and experiments with regards to facet formation 
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FIG. 1: (Color online). In-situ microscope images (a)-(g) of 
fatigue cracks in plexiglass at different stages of crack ad¬ 
vance in mixed mode I+III loading depicted schematically in 
(h) and corresponding example of crack-front segmentation 
in phase-field simulation (i). KiiijKi 0.3 in (a)-(e) and 

0.5 in (f)-(g); (a), (b) and (f) are experimental views from 
a direction approximately perpendicular to the plane of the 
parent crack with facets propagating downwards, while views 
(c), (d), (e) and (g) are views with the crack propagation di¬ 
rection out of the page. Views (c), (d) and (e) correspond 
to different stages of crack advance increasing from (c) to (e). 
Broken (pristine) regions of the samples appear in black (light 
blue) or darker (lighter) grey depending on the viewing direc¬ 
tion. The bar scale is 1 mm in all images. The red dashed 
lines in (a) highlight the curved fronts of two facets as guide 
to the eye; curved tips are clearly visible in (f). (i) Snapshots 
of phase-held fracture surfaces (0 = 1/2 surfaces) at different 
stages of crack advance increasing from top to bottom, show¬ 
ing that energetically favored A facets [18] propagate ahead 
of B facets eventually outgrowing them completely. Simula¬ 
tion parameters are GjGc — 1.5, KmlKi — 0.5, and box 
dimensions Dx = 307^, Dy = 100^ and Dz = 200^. 

and shed new light on the coarsening process. 

Experiments are carried out using plexiglas beams and 
a traditional three or four point bending setup [29]. To 
introduce some amount of mode III, the initial planar 
notch in the sample is tilted at an angle from the mode I 
central plane of symmetry [19, 30]. A special procedure 
is used to initiate a sharp crack with a straight front [29]. 
The corresponding values of the SIF for each angle and 
hence Km /Kj have been obtained by hnite element cal¬ 
culations, which show that KiiijKi varies between ap¬ 
proximately 0.1 and 0.5 when the notch angle varies be¬ 


tween 15° and 45°, where zero angle corresponds to pure 
mode I loading. This range was selected because it con¬ 
tains the linear instability threshold [Kiii/Ki)c ~ 0.39 
predicted by Eq. (2) for Poisson’s ratio of plexiglass 
V ~ 0.38. Finite element calculations also show [19] 
that KiiijKi is reasonably constant away from sample 
edges, thereby allowing us to investigate crack propaga¬ 
tion at constant KmjKi along a wide section of the 
parent crack inside the sample. Several beams were bro¬ 
ken by fatigue in the bending set-up [29]. The advantage 
of this cyclic type of loading is that the crack advance (i) 
is quasi-static, while leaving the crack path unchanged 
in comparison to the one obtained under monotonical in¬ 
creasing loading [31] and (ii) controlled by the number 
of cycles so that complex crack morphologies can be ob¬ 
served in-situ at different stages of crack growth. Obser¬ 
vations were made using a Leica binocular or a Keyence 
numerical microscope by transparency. 

Examples of experimental images are shown in Fig. 
l(a)-(g) for KmjKi values of 0.3 and 0.5 correspond¬ 
ing to initial notch angles of 30° and 45°, respectively. 
Those images reveal several important features. Firstly, 
facets have a finger-shape with curved tips and flat sides 
that is consistent with the shape predicted by phase- 
field simulations (Fig. l(i) and Movie 1 of [29]). Sec¬ 
ondly, facets form for values of Km /Kj both below and 
above the linear stability threshold {Km/Ki)c ~ 0.39. 
Within optical resolution, only energetically favored type 
A facets are observed to emerge from the parent crack 
with a well-defined tilt angle 0 from the original fracture 
plane. Thirdly, facets coarsen by elimination of other 
facets leading to an increase of both facet width and 
facet spacing along the array with increasing propagation 
length. Coarsening is clearly visible from top views in 
Fig. 1(b) and in the sequence Fig. 1(c)-(e), which more¬ 
over shows that surviving facets maintain the same angle 
while overgrowing others. Additional views are given in 
[29]. 

Simulations were carried out using a phase-field model 
of brittle fracture that, like gradient damage models 
[6, 32], regularizes stress-field divergences on a process 
zone scale ^ ^ around the crack front. All energy dissi¬ 
pation takes place on a characteristic timescale r [33]. As 
shown by an asymptotic analysis of the phase-field model 
in the limit where ^ is much smaller than all other dimen¬ 
sions [34], fracture in this model is governed by standard 
crack propagation laws assumed in the LEFM theoreti¬ 
cal framework, namely Griffith’s criterion and vanishing 
mode II SIF [35]. Since we are primarily interested in 
modeling crack evolution in a region away from the ex¬ 
perimental sample boundaries where Km /Kj is approx¬ 
imately uniform [19, 28], we carried out simulations in a 
rectangular slab geometry of length Dx^ width Dy and 
height Dz^ defined in Fig. 2(b), with the origin defined 
at the center of the slab. We impose fixed displacements 
at ^ = ±Dyl2^ Uy{x^±Dy/2^ z) = (mode I) and 
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FIG. 2: (Color online). Snapshots phase-field simulations 
illustrating the destabilization of planar crack growth for 
KiiijKi — 0.4. The crack propagation length a increases 
from (a) to (d) and both the crack front (blue lines) and its in¬ 
plane and out-of-plane projections (red lines) are shown, (e) 
Plot of linear instability threshold (Km jKi)^ versus DyjK. 
Planar growth is unstable (stable) above (below) the filled 
circles, where error bars reflect the uncertainty in stability 
threshold resulting from the fact that Km /Ki was increased 
in finite steps in the simulations. Km jKi values correspond¬ 
ing to the top (bottom) of each error bar were simulated and 
found to yield unstable (stable) propagation. In all simula¬ 
tions, G = 1.5Gc, = 230^ and = A = 60^. 

Uz{x^±Dy/2^ z) = (mode III), periodic boundary 

conditions in z that allow us to model a periodic array of 
daughter cracks infinite in 2; [24]. We use a “treadmill” 
that adds a strained (y^z) layer at x = Dx/2 and re¬ 
moves a layer at x = —D^ {2 when the crack has advanced 
by one lattice spacing. This allows us to simulate crack 
propagation lengths much longer than Dx (a ^ 
thereby modeling propagation in a slab infinitely long in 
X [29]. We also choose Dx > 2.bDy to eliminate the influ¬ 
ence of the two end-boundaries of the slab (x = :^Dxl2) 
on the central region of the slab (|x| ^ Dx) where the av¬ 
erage crack front position is maintained by the treadmill. 
Standard expressions of linear elasticity are used to re¬ 
late and A^ to the SIF [29] and therefore to KmlKj 
and GjGc where Gc ~ 27 (twice the surface energy) is 
known in the phase-field model [33, 34]. All simulations 
are performed with u = 0.38 of plexiglass. We simulated 
both quasi-static propagation, where the elastic field is 
relaxed at each time step of crack advance, and dynamic 
propagation by solving the full elastodynamic equations. 
Both sets of simulations yielded similar results for the 
range G/Gc < 1.5 where the ratio of the crack propa¬ 
gation speed to the shear wave speed vjc < 0.3 is small 
enough to neglect inertial effects [29]. 

We first carried out simulations to check quantita¬ 
tively the theoretical prediction of Eq. (2). For this 
purpose, we slightly perturbed the planar parent crack 
with a small amplitude helical perturbation of the form 
^^front 4“ ^^^front — ^0^ : whcrC ^^front S'lld ^^front 

indicate the x and y components of deviations of the 
front from the reference planar crack, respectively, and 
k = 2711Dz fits one wavelength Dz = oi the perturba¬ 
tion in the periodic domain in 2;. The stability of planar 


crack propagation is then determined by tracking the am¬ 
plitude of the perturbation that grows or decays exponen¬ 
tially in time [29] if propagation is unstable, as illustrated 
in Fig. 2(a)-(d), or stable, respectively. Simulations were 
carried out by increasing Km /Kj in small steps to deter¬ 
mine the threshold (Km jKi)c^ and repeating this pro¬ 
cedure for increasing values of Dy/K to quantify finite 
size effects. Fig. 2(e) shows that (Km/Ki)c increases 
monotonously with Dy/K and approaches a value rea¬ 
sonably close to the prediction (Km/Ki)c ~ 0.39 of Eq. 
(2) in the large system size (Dy/K 1) limit. Con¬ 
sistent with the result of Eig. 2(e), an examination of 
strain fields shows that finite size effects becomes neg¬ 
ligible when Dy/K > 2 [29]. We conclude that LEEM 
theory (Eq. (2)) and phase-field modeling predict simi¬ 
lar linear instability thresholds in the large system size 
limit, and that facets are experimentally observed well 
below this threshold. 

Next, in order to explore the nonlinear character of 
the bifurcation from planar to segmented crack front, we 
measured experimentally the facet tilt angle d extracted 
from three-dimensional maps of post-mortem fracture 
surfaces obtained using a profilometer as detailed in [28]. 
The angle d is plotted versus Km/Kj in Eig. 3(a). 
Eurthermore, we investigated computationally the prop¬ 
agation of periodic arrays of A facets, where coarsen¬ 
ing is suppressed by choosing Dz = K due to the peri¬ 
odic boundary conditions along z. In this geometry, we 
tracked the steady-state branch of propagating solutions 
by decreasing Km/Kj starting from values above the 
linear instability threshold to values below this thresh¬ 
old, as low as 0.07 to span the entire experimental range 
of mode mixity. Eor each Km /Ki value, we allowed the 
facet to relax to a new stationary shape and tilt angle, as 
illustrated in Eig. 3(b) for a simulation where Km/Kj 
was decreased from 0.5 to 0.07. The computed tilt angles 
are compared to experimental results in Eig. 3(a) with 
the corresponding facet shapes shown in Eig. 3(c). Both 
the facet shapes, which gently curve at their extremi¬ 
ties in the yz plane due to elastic interactions between 
neighboring facets, and the tilt angles are in good quanti¬ 
tative agreement with experimental observations within 
measurement errors. Eig. 3(a) also shows that computed 
tilt angles are weakly dependent on system size (Dy/K) 
and fall below the prediction of a simple theory, which 
assumes that facets are shear-free [16, 24]. Those results 
demonstrate that propagating segmented front solutions 
exist over the entire range of Km/Kj investigated ex¬ 
perimentally, including values less than (Km/Ki)^. We 
conclude that the bifurcation from planar to segmented 
front is strongly subcritical, with bistability of planar and 
segmented crack growth for Km/Kj < (Km/Ki)c as 
illustrated schematically in Eig. 3(d). 

To characterize coarsening in phase-field simulations, 
we investigated the stability of periodic array of facets by 
repeating the above series of simulations with two facets 
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FIG. 3: (Color online), (a) Comparison of facet tilt an¬ 
gles obtained from experiments and simulations, where red 
and blue arrows indicate the instability thresholds of planar 
crack propagation for Dy/A = 1 and Dy/A = 2, respectively 
(see Fig. 2(e)), and theoretically predicted assuming shear- 
free facets (dashed line) [16, 24]. (b) Snapshots of a phase- 
field simulation for Dy/A — 1 demonstrating the subcriti- 
cal nature of the bifurcation from planar to segmented crack 
propagation. A propagating segmented front solution for 
Kiii/Ki = 0.5 was used as initial condition (Q — 31°). The 
facet continuously rotated towards a lower angle in response 
to the decrease in Km /Ki and then reached its steady state 
(Q — 11.2°) after propagating a distance a — 5.5A (see Movie 
2 of [29]). (c) Out-of-plane and in-plane (inset) crack-front 
projections. In all simulations, D^ — 154^, Dy = Dz — 60^, 
A = 60^ and G — 1.5Cc. (d) Schematic diagram of sub- 

critical bifurcation recapitulating the experimental and sim¬ 
ulations results with solid (dashed) lines representing stable 
(unstable) solutions. 

{Dz = 2A). This geometry is motivated by the striking 
similarity between the coarsening behavior of facets in 
the present experiments (Fig. l(a)-(g)) and coarsening 
of curved fronts in other interfacial pattern forming sys¬ 
tems, in particular viscous fingering [36] and dendritic 
crystal growth [37, 38]. In those systems, it is well- 
established that coarsening of finger arrays is associated 
with a spatial period-doubling linear instability of the ar¬ 
ray leading to elimination of one of every two fingers in 
the array by exponential amplification of small perturba¬ 
tions. Results of simulations illustrated in Fig. 4(a) show 
that arrays of facets exhibit a similar spatial period dou¬ 
bling instability driven by elastic interactions between 



z/A 


FIG. 4: (Color online), (a) Illustration of spatial period dou¬ 
bling instability in a phase-field simulation for Km/Kj = 
0.5; out-of-plane and in-plane projections of crack fronts at 
different times are plotted in the top panel and the bottom 
panel, respectively (see Movie 3 of [29]). (b) Semi-log plot 
of difference of tip positions along the propagation a:-axis be¬ 
tween leading and lagging facets versus scaled time for dif¬ 
ferent Km/Kj. Inset: coarsening rate /3 versus Km/Kj 
obtained from experiments and phase-field simulations. In all 
simulations, Dx — 307^, Dy = 60^, Dz = 120^, A = 60^ and 
C = 1.5Cc. 

facets. This instability yields an increase (decrease) of 
the SIF and hence the energy release rate at the tips of 
leading (lagging) facets. The amplification rate of insta¬ 
bility is obtained by computing the difference of x-tip po¬ 
sition Axtip(t) between leading and lagging facets, which 
grows exponentially in time starting from an infinitesi¬ 
mal perturbation, Axtip(t) Axtip(0)e^'^o^/^, where vq 
and A are the initial facet growth velocity and spacing, 
respectively. The slopes of semi-log plots of Axtip(t)/A 
versus v^t/A in Fig. 4(b) yield values of uj that increase 
markedly with Km/Kj^ showing that a larger mode III 
component leads to a faster elimination rate of facets. 

Coarsening, clearly visible in Fig. 1(b) and other ex¬ 
perimental views [29], was quantified experimentally by 
analyzing post-mortem fracture surfaces [28]. The results 
show that the relation between the mean facet spacing 
A and the crack propagation length a is approximately 
linear, with a mean slope P = dA/da increasing with 
KiiilKi (inset of Fig. 4(b)). To relate the coarsen¬ 
ing rates in phase-field simulations and experiments, we 
derive a simple evolution equation for the average ar¬ 
ray spacing A based on dynamical mean-field picture as 
previously done for dendritic arrays [37]. The coarsen¬ 
ing rate p = dA/da ~ AA/Aa where A A is the change 
of array spacing due to elimination of one of every two 
facets along the array or AA « A, while Aa is the dis¬ 
tance that the facets propagated during the elimination 
process. Since elimination occurs via exponential ampli¬ 
fication of small perturbations, facets will propagate an 
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average distance Aa ^ K/uj during this process, yielding 
the prediction ^ cj, or = Cuo where C is a constant 
prefactor of order unity. The comparison in the inset of 
Fig. 4(b) shows that this simple theory is able to pre¬ 
dict reasonably well the increase of the coarsening rate 
with KiiijKi up to the value of the constant prefac¬ 
tor C — 0.198 determined from a global best fit to the 
experimental data for all KnijKi values. 

The present results reconcile the prediction of linear 
stability analysis (Eq. (2)) with experimental observa¬ 
tions by showing that the bifurcation from planar to 
segmented crack growth is strongly sub critical; facet ar¬ 
rays exist as fundamental crack propagating solutions of 
LEFM for a range of KmjKi values extending below 
the instability threshold. They further show that coars¬ 
ening is driven by a spatial period doubling instability of 
facet arrays with a growth rate that depends on mode 
mixity. The reasonably good quantitative agreement be¬ 
tween simulated and observed morphologies suggests that 
LEFM is an adequate theory to describe complex geomet¬ 
rical features of both brittle and fatigue cracks in mixed 
mode I+III fracture. While the present results show that 
the subcritical propagation of segmented cracks is theo¬ 
retically possible, they do not identify the mechanism 
and scale of subcritical facet formation. As suggested 
by a recent LEFM analysis, materials imperfections may 
contribute to this process [39]. However, this scenario, 
and even more fundamentally the ability of LEFM to 
model subcritical facet formation, remain to be explored 
both computationally and experimentally. 
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EXPERIMENTS 



FIG. 1: Left hand side: 3 Point Bending setup. Right hand side: Lateral view of the initial slit made by micro-milling and of 
the additional sharp crack introduced by pushing under a controlled manner, a razor blade in the initial slit. 


Experiments were carried out using plexiglas beams of dimensions L x W x b (Fig. 1, left). We used two different 
sizes corresponding to L = 260 mm, LF = 60 mm, 6 = 10 mm for the large ones and L = 100 mm, LF = 10 mm, 
b = 10 mm for the small ones. The initial slit with a blunted tip of radius R = 300 /am was made by micro-milling. 
To initiate a true crack with a sharp tip and a smooth straight front (Fig. 1, right), we pushed a wedge (razor blade) 
into the slit, quasi-statically, under controlled slowly increasing force, using a tensile machine. The total length d of 
the crack (slit+sharp crack) is d ^ LF/3 in both cases. The residual stresses introduced by the slit manufacturing, 
were relaxed by heating the samples at 90-95 °C during 10 hours. Their annihilation is checked using polarizers. To 
introduce some amount of mode III, the initial planar notch in the sample is tilted at an angle Fq from the mode I 
central plane of symmetry, Fq varying between 15° and 45°, where zero angle corresponds to pure mode I loading. 
Larger and smaller samples were loaded in 3 and 4 point bending setups, respectively. The loading frequency was 
/ = 5 Hz and Kminf Kmax = 0-1. The amplitude of the stress intensity factor AK ^ 0.5 MPa-m^/^ has been chosen 
well below the brittle fracture threshold Kc ^ 1 MPa-m^/^, to avoid brittle fracture, while being large enough to 
ensure propagation [1]. 

Sets of experiments were performed on series of similar beams by stopping the propagation at different stages of the 
propagation. Figure 2 gives some in-situ views of the corresponding crack morphologies, acquired with a numerical 
Keyence microscope. Front and bottom views are given in the first and second row respectively. Each column 
corresponds to the same sample. The beam size is L = 100 mm, LF = 10 mm, b = 10 mm; columns (a)-(c) correspond 
to Fq = 30° {Kiii/Ki « 0.3) and column (d) to Fq = 20° {Kjji/Kj « 0.18). The development of rotated facets that 
coalesce during propagation are clearly visible. One can observe i) comparing columns (a) to (c), that the rotation 
angle of those facets is approximately constant during propagation, for a given value of KiiijKi and ii) comparing 
columns (a)-(c) with (d), that the rotation angle of these facets decreases with decreasing KinjKi. 

The coalescence rate P cannot be quantified on the in-situ samples due to optical distortions induced by the plexiglas. 
For this purpose, other samples were broken completely (Fig. 3 (a)-(c)) and three dimensional profilometer maps of 
their fracture surfaces were done. Quantification of 0 (Fig. 3(a) of the main text) and p (inset of Fig. 4(b) of the main 
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FIG. 2: In-situ numerical microscope images of partially broken samples. Each column corresponds to one sample, the first 
row being a front view and the second row being a bottom view as indicated by cartesians axes corresponding to Fig. 1(h) of 
the main text. KmlKi ^ 0.3 in (a)-(c) and ~ 0.18 in (d). In the front views, the initial slit and the facets appear in black. 
In the bottom views, the initial larger slit is dark and the facets appear in white. The bar scale is 1 mm in all images. 
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FIG. 3: Post-mortem front views of the fracture surfaces corresponding to (a) KiiijKi — 0.15, (b) KiiijKi — 0.3, (c) 
KiiilKi — 0.5 (L = 260 mm, VF = 60 mm, 5 = 10 mm). The white zones corresponds to the partially broken fragments left 
behind on the crack surface. The formation of these fragments is due to the lateral propagation of two adjacent interacting 
facets as sketched on (d). The bar scale is 2 mm in (a)-(c). 

text) has been done performing some statistical post-treatment over several facets of these profiles. The procedure is 
explained in detail in [2]. 

The more clear zones in Fig. 3(a)-(c) correspond to fragments left behind by partial breaking of the zones located 
between adjacent facets. Indeed, beside propagating in the x—direction, the facets propagate also in the lateral 
z—direction. As they interact, they curve to form the well-known [3, 4] “en-passant S-shape” pattern (see Fig. 3(d)). 
Since the symmetry of this interaction is sensible to any perturbation, generally only one out of two crack tips connects 
to the adjacent facet, leaving on one of the fracture surface, a partially broken fragment. These fragments are visible 
in white on Fig. 3(a)-(c) and are lacking on the complementary fracture surface. 


SIMULATIONS 

We conducted large-scale simulations of mixed-mode I+III fracture using the phase-field model originally proposed 
by Karma, Kessler and Levine (KKL) [5], which has been used to solve fracture propagation problems over the last 
decade [6-10]. The model introduces a scalar order parameter 0 to distinguish between intact and broken states of 
the material. The total energy of the system is given by the functional [5, 11] 

^ = /{f(^) + 9 (0) (estrain - 60) | dV, ( 1 ) 

where p is the mass density; g (0) = 40^ — 30^ is conventionally chosen with the properties that ^ (0) = 0, ^ (1) = 1 
and g' (1) = g' (0) = 0 [5, 11]; u = {ux^Uy^Uz) represents the displacement field; estrain = A (e^i)^ /2 + /i is the 
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strain energy density and Sij = {diUj djUi) /2 is the strain tensor of a linear elastic material with i = 1, 2 and 3 
corresponding to y and z, respectively, A and y are the Lame coefficients. The function g decreases monotonously 
from the intact state (0 = 1) to the fully broken state (0 = 0) and accounts for elastic softening at large strain. When 
the strain energy density egtrain exceeds the threshold Cc, the broken state becomes energetically favored. Energy 
dissipation occurs in the process zone of size ^ tz/ (2ec) around the crack tip, where 0 varies smoothly between 
0 and 1. In the phase-field model, the Griffith’s criterion Gc is given by Gc = 0^1 — g (0) [5, 11]; the 

iso-surfaces of cj) = 0.5 are conventionally defined as the fracture surfaces. The equations of motion for the phase-field 
(j) and displacement field u are derived variationally from the energy functional [5, 11]: 
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We studied crack-front instabilities of mixed mode I+III fracture propagation along the x-axis in a rectangular slab 
of length Dx, width Dy and thickness Dz with center at the origin. The mixed I+III loading was imposed by fixed 
displacements at ^ = EDyj2^ Uy (x, EDyj^^ z) = (mode I) and Uz {x, EDyl2, z) = (mode III), and periodic 
boundary conditions in the ^-direction. To allow full relaxation of the crack front to a stationary state, simulations 
were carried out on a “treadmill” that adds a strained {y^z) layer at x = Dx/2 and removes a layer at x = —0^12 
when the crack has advanced by one lattice spacing. This allows us to simulate crack propagation lengths much longer 
than Dx {a ^ Dx), thereby modeling propagation in a slab infinitely long in x. 

The crack dynamics is controlled by two key parameters KjjjjKj and GjGc, where Kj = 2fi/S.y^j2 (1 + A///) /Dy 
and Km = 2iil\z/^J^y are the stress intensity factors of a semi-infinite planar crack and G = 
((1 — u) Kj + Kjjj) / ( 2 //) is the corresponding energy release rate, where z/ is Poisson’s ratio. Simulations were 
carried out with the range G/Gc < 1.5 where the ratio of the crack propagation speed to the shear wave speed 
v/c < 0.3, which is small enough to neglect inertial effects. The slab length Dx was chosen greater than 2.5 of the 
slab width Dy to eliminate the influence of the two end-boundaries of the slab {x = EDxl2) on the central region of 
the slab {\x\ Dx) where the average crack front position is maintained by the treadmill. As a result, the planar 
crack propagation speed becomes independent on the system length, as shown in Fig. 4. 

We rewrite the phase-field model in a dimensionless form by measuring length in units of the fracture process zone 
scale ^ and time in units of the characteristic dissipation time scale r = 1/ (/ry). The remaining simulation parameters 
are chosen as follows, Cc/g = 1 / 2 , u = 0.38 and cr/^ = 2.76, where the shear wave speed c = g/p. Large scale 
simulations of the order of 10^ — 10^ grid points were performed using graphics processing units (GPUs) with the 
GUDA parallel programming language. 

The energy functional of the phase-field model was discretized in space on a cubic grid of uniform mesh size A = 0.3^ 
(Fig. 5). To write the discretized energy in compact form, we define the superscripts m, and n to refer to the 
gridpoint at x = ^A, y = mA, and z = nA, the subscripts = {1,2,3} to refer to the corresponding {x^y^z} 

components of the displacement fields, and use the standard Kronecker delta defined by 5ij = 1 \ii = j and 5ij = 0 if 
i 7 ^ j. We discretize duijdxj^ {duijdx^^ and (90/5xj)^ on the gridpoint by the following approximations: 
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Accordingly, the strain energy density egtrain on gridpoint {£, m, n) becomes 
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FIG. 4: Plot of pure mode I crack propagation speed v (scaled by the shear wave speed c) versus GjGc- Phase-held simulations 
show that the crack speed becomes independent of the simulation box length when DxjDy > 2.5 (inset). Simulation parameters 
are Dy = 60^ and Dz = 30^. 
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The discretized form of the energy functional Eq. (1) on the cubic grid becomes 
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where, in addition to the functions dehned above, 

p£, 7 n, 7 i _ ^^£+6ii,7n+62i,7i+63i _ ^£, 77 i, 7 i^‘^ _j_ ^^£, 7 n, 7 i _ ^£— 6 ii, 7 n—621,71—63i^‘^ 

From Eq. (2), Eq. (3) and Eq. (4), we obtain the equations of motion 
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FIG. 5: Spatial discretization of phase-field model on a cubic grid of uniform mesh size A. Gridpoint labeled 
represents a unit cubic of material of mass pA^ at location (^A,mA,nA). 



FIG. 6: (Golor online). Semi-log plots of perturbation amplitude A versus scaled time in simulations. The perturbation 
amplitude A exponentially grew in time for KmlKi = 0.35 (red squares) and decayed in time for KmlKi = 0.25 (blue 
circles). Simulation parameters are G = 1.5Gc, Aq = 0.6^, Dx = 230^, Dy = 120^ and Dz = 30^. The same linear instability 
threshold was obtained in simulations using Aq = 1.2^ and Aq = 1.8^. 

^ a(E/A°) 

which are integrated in time using a Verlet scheme with a timestep size At/r = 0.001. 

To investigate the onset of instability of mixed mode I+III crack propagation, we carried out simulations starting 
from a parent planar crack in yz plane for x <0. As in [10], the planar crack front was initially perturbed by a helical 
perturbation, Sx^ront + ^^^front = where Sx^ront and ^^uont indicate the x and y components of deviations 

of the front from the reference planar crack, respectively. The amplitude A (t) = was then found to grow 

(<j > 0) or decay (a < 0) exponentially in time (Fig. 6). 

The theoretical instability threshold predicted by Eq. (2) of the main text is only valid in the large system size limit 
[12]. To characterize finite size effects, simulations with DyjA between 0.5 and 4 were carried out. The results show 
that finite size effects reduce significantly [KiiijKi)^ when Dy is comparable to A (Fig. 2(e) of the main text). For 
such system sizes, the strain fields along the y = -ADyj2 boundaries are spatially varying along 2 ; and hence strongly 
influenced by the facets. However for Dy/K > 2, the strain fields become nearly independent of 2 ; as expected in the 
large system size limit Dy ^ A (Fig. 7). In addition, the results of simulations with A/^ between 30 and 60 show 
that the crack-front instability is independent of A in the limit A 
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FIG. 7: Finite size effects: (a) (Color online). Images of strain field Syz of a slice along x behind the facet tips (the corresponding 
slice position is illustrated in the top panel) from simulations with = 1, 2 and 4. (b) Plots of strain field Syz of the sliced 

regions oX y — Dyjl versus z (scaled by A). The strain fields on boundaries are strongly infiuenced by facets when DyjK — 1 
(red line) and nearly constant along z for both DyjK — 2 (blue dashed line) and DyjA. — 4 (green double dashed line), where 
boundary effects become negligible. In all simulations, KmlKi — 0.35, DxjDy = 3.75 and Dz = A = 30^. 
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